#### setting environment ####
require(quanteda)
require(stm)
require(runjags)
require(coda)

newspaper.names <- c("Asahi", "Chugoku", "Chunichi", "Hokkaido", "Kahoku", 
                     "Mainichi", "Nikkei", "Nishinippon", "Sankei", "Yomiuri")

#### data preparation ####
# estimation results of the Wordfish models
load("Study2_Wordfish_result.Rdata")

#### preliminary factor analysis to impose a sign restriction ####
# function to create initial values
initial.fun.preliminary <- function(seed1, seed2) {
  set.seed(seed1)
  b <- rnorm(65)
  x <- rnorm(10, 0, 0.01)
  p <- runif(65, 0.5, 1)
  list(beta = b, x = x, inv.phi = p, 
       .RNG.name = "base::Wichmann-Hill", .RNG.seed = seed2)
}
# create initial values
initial.values.preliminary <- initial.fun.preliminary(12345, 67890)
# estimate a factor analysis model using JAGS
FA.result.preliminary <- run.jags(model = "factor_analysis_JAGS_preliminary.R", 
                                    monitor = c("x", "beta", "phi"), 
                                    data = list(y = topic.specific.positions, n = 10, m = 65), 
                                    inits = initial.values.preliminary, 
                                    n.chains = 1, burnin = 5000, sample = 1000, modules = "glm", 
                                    summarise = FALSE, thin = 10, method = "parallel")

## fix the signs of factor loadings of some items so that Chunichi's ideal point becomes negative
# whether Chunichi's ideal point was estimated to be negative
Chunichi.score.sign.preliminary <- sign(mean(as.matrix(FA.result.preliminary$mcmc[, 3]))) == -1
# posterior means of factor loadings
post.loadings.preliminary <- colMeans(as.matrix(FA.result.preliminary$mcmc[, 11:75]))
# IDs of five topics whose factor loadings should be positive to make Chunichi's ideal point negative
positive.items <- order(post.loadings.preliminary, decreasing = Chunichi.score.sign.preliminary)[1:5]
# IDs of other topics
unconstrained.items <- c(1:65)[-positive.items]

#### main factor analysis ####
# function to create initial values
initial.fun <- function(seed1, seed2) {
  set.seed(seed1)
  b <- rnorm(65)
  b[positive.items] <- runif(5, 10, 15)
  x <- rnorm(10, 0, 0.01)
  p <- runif(65, 0.5, 1)
  list(beta = b, x = x, inv.phi = p,
       .RNG.name = "base::Wichmann-Hill", .RNG.seed = seed2)
}
# create initial values
initial.values <- list()
for (i in 1:3) {
  initial.values[[i]] <- initial.fun(100 * i, 1000 * i)
}
# estimate a factor analysis model using JAGS
FA.result <- run.jags(model = "factor_analysis_JAGS.R", 
                        monitor = c("x", "beta", "phi"), 
                        data = list(y = topic.specific.positions, n = 10, m = 65, 
                                    positive.items = positive.items, 
                                    unconstrained.items = unconstrained.items), 
                        inits = initial.values, 
                        n.chains = 3, burnin = 5000, sample = 5000, modules = "glm", 
                        summarise = FALSE, thin = 20, method = "parallel")
# convergence check
which(gelman.diag(FA.result$mcmc, multivariate = FALSE)[[1]][, 1] > 1.1)

# save(FA.result, file = "FA_result.Rdata")